home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / gammabase.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  7KB  |  303 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. /* forward declarationss */
  15. #ifdef ANSI
  16. double gammp(double,double),gnorm(double,double),gser(double,double,double),
  17.        gcf(double,double,double),gammad(double *,double *,int *),
  18.        ppchi2(double *,double *,double *,int *);
  19. #else
  20. double gammp(),gnorm(),gser(),gcf),gammad(),ppchi2();
  21. #endif ANSI
  22.  
  23. void gammabase(x, a, p)
  24.      double *x, *a, *p;
  25. {
  26.   *p = gammp(*a, *x);
  27. }
  28.  
  29. double ppgamma(p, a, ifault)
  30.      double p, a;
  31.      int *ifault;
  32. {
  33.   double x, v, g;
  34.  
  35.   v = 2.0 * a;
  36.   g = gamma(a);
  37.   x = ppchi2(&p, &v, &g, ifault);
  38.   return(x / 2.0);
  39. }
  40.  
  41. /*
  42.   Static Routines
  43. */
  44.  
  45. /*
  46.   From Numerical Recipes, with normal approximation from Appl. Stat. 239
  47. */
  48.  
  49. #define EPSILON 1.0e-14
  50. #define LARGE_A 10000.0
  51. #define ITMAX 1000
  52.  
  53. static double gammp(a, x)
  54.      double a, x;
  55. {
  56.   double gln, p;
  57.  
  58.   if (x <= 0.0 || a <= 0.0) p = 0.0;
  59.   else if (a > LARGE_A) p = gnorm(a, x);
  60.   else {
  61.     gln = gamma(a);
  62.     if (x < (a + 1.0)) p = gser(a, x, gln);
  63.     else p = 1.0 - gcf(a, x, gln);
  64.   }
  65.   return(p);
  66. }
  67.  
  68. /* compute gamma cdf by a normal approximation */
  69. static double gnorm(a, x)
  70.      double a, x;
  71. {
  72.   double p, sx;
  73.  
  74.   if (x <= 0.0 || a <= 0.0) p = 0.0;
  75.   else {
  76.     sx = sqrt(a) * 3.0 * (pow(x / a, 1.0 / 3.0) + 1.0 / (a * 9.0) - 1.0);
  77.     normbase(&sx, &p);
  78.   }
  79.   return(p);
  80. }  
  81.   
  82. /* compute gamma cdf by its series representation */
  83. static double gser(a, x, gln)
  84.      double a, x, gln;
  85. {
  86.   double p, sum, del, ap;
  87.   int n, done = FALSE;
  88.  
  89.   if (x <= 0.0 || a <= 0.0) p = 0.0;
  90.   else {
  91.     ap = a;
  92.     del = sum = 1.0 / a;
  93.     for (n = 1; ! done && n < ITMAX; n++) {
  94.       ap += 1.0;
  95.       del *= x / ap;
  96.       sum += del;
  97.       if (fabs(del) < EPSILON) done = TRUE;
  98.     }
  99.     p = sum * exp(- x + a * log(x) - gln);
  100.   }
  101.   return(p);
  102. }
  103.  
  104. /* compute complementary gamma cdf by its continued fraction expansion */
  105. static double gcf(a, x, gln)
  106.      double a, x, gln;
  107. {
  108.   double gold = 0.0, g, fac = 1.0, b1 = 1.0;
  109.   double b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
  110.   double p;
  111.   int done = FALSE;
  112.  
  113.   a1 = x;
  114.   p = 0.0;
  115.   for(an = 1.0; ! done && an <= ITMAX; an += 1.0) {
  116.     ana = an - a;
  117.     a0 = (a1 + a0 * ana) * fac;
  118.     b0 = (b1 + b0 * ana) * fac;
  119.     anf = an * fac;
  120.     a1 = x * a0 + anf * a1;
  121.     b1 = x * b0 + anf * b1;
  122.     if (a1 != 0.0) {
  123.       fac = 1.0 / a1;
  124.       g = b1 * fac;
  125.       if (fabs((g - gold) / g) < EPSILON) {
  126.         p = exp(-x + a * log(x) - gln) * g;
  127.         done = TRUE;
  128.       }
  129.       gold = g;
  130.     }
  131.   }
  132.   return(p);
  133. }
  134.  
  135. static double gammad(x, a, iflag)
  136.      double *x, *a;
  137.      int *iflag;
  138. {
  139.   double cdf;
  140.  
  141.   gammabase(x, a, &cdf);
  142.   return(cdf);
  143. }
  144.  
  145. /*
  146.   ppchi2.f -- translated by f2c and modified
  147.  
  148.   Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35
  149.   To evaluate the percentage points of the chi-squared
  150.   probability distribution function.
  151.   
  152.   p must lie in the range 0.000002 to 0.999998,
  153.     (but I am using it for 0 < p < 1 - seems to work)
  154.   v must be positive,
  155.   g must be supplied and should be equal to ln(gamma(v/2.0)) 
  156.   
  157.   Auxiliary routines required: ppnd = AS 111 (or AS 241) and gammad.
  158. */
  159.  
  160. static double ppchi2(p, v, g, ifault)
  161.      double *p, *v, *g;
  162.      int *ifault;
  163. {
  164.   /* Initialized data */
  165.  
  166.   static double aa = .6931471806;
  167.   static double six = 6.;
  168.   static double c1 = .01;
  169.   static double c2 = .222222;
  170.   static double c3 = .32;
  171.   static double c4 = .4;
  172.   static double c5 = 1.24;
  173.   static double c6 = 2.2;
  174.   static double c7 = 4.67;
  175.   static double c8 = 6.66;
  176.   static double c9 = 6.73;
  177.   static double e = 5e-7;
  178.   static double c10 = 13.32;
  179.   static double c11 = 60.;
  180.   static double c12 = 70.;
  181.   static double c13 = 84.;
  182.   static double c14 = 105.;
  183.   static double c15 = 120.;
  184.   static double c16 = 127.;
  185.   static double c17 = 140.;
  186.   static double c18 = 1175.;
  187.   static double c19 = 210.;
  188.   static double c20 = 252.;
  189.   static double c21 = 2264.;
  190.   static double c22 = 294.;
  191.   static double c23 = 346.;
  192.   static double c24 = 420.;
  193.   static double c25 = 462.;
  194.   static double c26 = 606.;
  195.   static double c27 = 672.;
  196.   static double c28 = 707.;
  197.   static double c29 = 735.;
  198.   static double c30 = 889.;
  199.   static double c31 = 932.;
  200.   static double c32 = 966.;
  201.   static double c33 = 1141.;
  202.   static double c34 = 1182.;
  203.   static double c35 = 1278.;
  204.   static double c36 = 1740.;
  205.   static double c37 = 2520.;
  206.   static double c38 = 5040.;
  207.   static double zero = 0.;
  208.   static double half = .5;
  209.   static double one = 1.;
  210.   static double two = 2.;
  211.   static double three = 3.;
  212.  
  213. /*
  214.   static double pmin = 2e-6;
  215.   static double pmax = .999998;
  216. */
  217.   static double pmin = 0.0;
  218.   static double pmax = 1.0;
  219.  
  220.   /* System generated locals */
  221.   double ret_val, d_1, d_2;
  222.   
  223.   /* Local variables */
  224.   static double a, b, c, q, t, x, p1, p2, s1, s2, s3, s4, s5, s6, ch;
  225.   static double xx;
  226.   static int if1;
  227.  
  228.  
  229.   /* test arguments and initialise */
  230.   ret_val = -one;
  231.   *ifault = 1;
  232.   if (*p <= pmin || *p >= pmax) return ret_val;
  233.   *ifault = 2;
  234.   if (*v <= zero) return ret_val;
  235.   *ifault = 0;
  236.   xx = half * *v;
  237.   c = xx - one;
  238.  
  239.   if (*v < -c5 * log(*p)) {
  240.     /* starting approximation for small chi-squared */
  241.     ch = pow(*p * xx * exp(*g + xx * aa), one / xx);
  242.     if (ch < e) {
  243.       ret_val = ch;
  244.       return ret_val;
  245.     }
  246.   }
  247.   else if (*v > c3) {
  248.     /* call to algorithm AS 111 - note that p has been tested above. */
  249.     /* AS 241 could be used as an alternative. */
  250.     x = ppnd(*p, &if1);
  251.  
  252.     /* starting approximation using Wilson and Hilferty estimate */
  253.     p1 = c2 / *v;
  254.     /* Computing 3rd power */
  255.     d_1 = x * sqrt(p1) + one - p1;
  256.     ch = *v * (d_1 * d_1 * d_1);
  257.  
  258.     /* starting approximation for p tending to 1 */
  259.     if (ch > c6 * *v + six)
  260.       ch = -two * (log(one - *p) - c * log(half * ch) + *g);
  261.   }
  262.   else{
  263.     /* starting approximation for v less than or equal to 0.32 */
  264.     ch = c4;
  265.     a = log(one - *p);
  266.     do {
  267.       q = ch;
  268.       p1 = one + ch * (c7 + ch);
  269.       p2 = ch * (c9 + ch * (c8 + ch));
  270.       d_1 = -half + (c7 + two * ch) / p1;
  271.       d_2 = (c9 + ch * (c10 + three * ch)) / p2;
  272.       t = d_1 - d_2;
  273.       ch -= (one - exp(a + *g + half * ch + c * aa) * p2 / p1) / t;
  274.     } while (fabs(q / ch - one) > c1);
  275.   }
  276.  
  277.   do {
  278.     /* call to gammad and calculation of seven term Taylor series */
  279.     q = ch;
  280.     p1 = half * ch;
  281.     p2 = *p - gammad(&p1, &xx, &if1);
  282.     if (if1 != 0) {
  283.       *ifault = 3;
  284.       return ret_val;
  285.     }
  286.     t = p2 * exp(xx * aa + *g + p1 - c * log(ch));
  287.     b = t / ch;
  288.     a = half * t - b * c;
  289.     s1 = (c19 + a * (c17 + a * (c14 + a * (c13 + a * (c12 + c11 * a))))) / c24;
  290.     s2 = (c24 + a * (c29 + a * (c32 + a * (c33 + c35 * a)))) / c37;
  291.     s3 = (c19 + a * (c25 + a * (c28 + c31 * a))) / c37;
  292.     s4 = (c20 + a * (c27 + c34 * a) + c * (c22 + a * (c30 + c36 * a))) / c38;
  293.     s5 = (c13 + c21 * a + c * (c18 + c26 * a)) / c37;
  294.     s6 = (c15 + c * (c23 + c16 * c)) / c38;
  295.     d_1 = (s3 - b * (s4 - b * (s5 - b * s6)));
  296.     d_1 = (s1 - b * (s2 - b * d_1));
  297.     ch += t * (one + half * t * s1 - b * c * d_1);
  298.   } while (fabs(q / ch - one) > e);
  299.  
  300.   ret_val = ch;
  301.   return ret_val;
  302. }
  303.